O intuito desta páginas é demonstrar as avaliações feitas para validação do modelo elaborado para a dissertação de mestrado no programa PROCAM-USP
O modelo hidrodinâmico é a base que fornecerá ao seguimento do estudo, a possibilidade de entender e simular os processos costeiros e como estes impactam nos compartimentos bióticos, sociais e econômicos
Para isso, foi considerada extensão desde o Município de São Sebastião Até o limite do Munícipio de Ubatuba. Nas figuras a seguir são demonstradas a grade e a extensão do domínio numérico.
A grade contem elementos uniformemente distribuidos em dimensões horizontais. Cada célula é de aproximadamente 50 x 50 m de largura.
Para a discretização vertical, foram utilizadas 8 camadas SIGMA, visto que na região, conforme a literatura, há presença marcada de termoclina no verão devido à intrusão da ACAS. Além disso, há variação nas profundidades de 60 - 0 m, sendo assim, o sistema de coordenadas sigma verticais é o mais adequado para representar essa variação no domínio.
O Modelo utilziado foi o Delft3D em seu modulo hidrodinâmico.
Os parametros considerados foram:
Foi utilizada a biblioteca Delft Dashboard para gerar as fronteiras abertas e as constantes harmônicas do TPXO 8.0. Segue imagem das fronteiras:
Sendo assim, 6 fronteiras a Oeste, 5 fronteiras a Sul e 2 fronteiras a Leste.
Em cada um dos domínios abertos, além das constantes harmonicas, foram adicionados os dados disponíveis no sistemas de modelo HYCOM. Pela disponibilidade de nível e corrente, foram utilizadas bordas Riemann, que são bordas não reflexivas, que consideram propagam snal de velocidade, levando em conta as variações do nível e a batimetria local:
Cada borda então conta com:
Devido a extensão do domínio, ao longo dos testes, verificou-se que apenas a inclusão de forçantes oceanográficas não era o suficiente. Para complementar, foi acionada a troca de calor no modelo, utilzando as equações referentes ao modo "ocean". Para implementação deste sistema, foi preciso:
Felizmente, todas as variáveis descritas estão disponíveis no sistema de reanálises CFSR. Desta maneira foram utilizadas tais componentes, além do vento e pressão. Todos mapeados ao longo do domínio do modelo.
As forçantes, por serem de outros modelos, não são uma restrição para escolha do período de modelagem, no entanto a ausencia de dados sim. Segue imagem de todas as informações fornecidas pelo BNDO:
Mesmo com uma infinidade de dados do BNDO, apenas um conjunto foi longevo suficiente para análise. Os dados utilizados, são referentes a uma coleta de janeiro de 2014 até final de maio de 2014, na baía do Araçá.
O conjunto de dados apresenta:
Assim, como a análise de constantes harmônicas não é um limitante, decidiu-se para calibração rodar de: 01/01/2014 até 01/06/2014. Compreendendo assim parte do verão e outono.
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import xarray as xr
from datetime import timedelta
from TTutils.Delft3D import *
from TTutils.MO import *
from TTutils import Calibration
import matplotlib.pyplot as plt
import matplotlib.dates as mp_dates
SMALL_SIZE = 8
MEDIUM_SIZE = 10
BIGGER_SIZE = 12
plt.rc('font', size=BIGGER_SIZE) # controls default text sizes
plt.rc('axes', titlesize=BIGGER_SIZE) # fontsize of the axes title
plt.rc('axes', labelsize=BIGGER_SIZE) # fontsize of the x and y labels
plt.rc('xtick', labelsize=BIGGER_SIZE) # fontsize of the tick labels
plt.rc('ytick', labelsize=BIGGER_SIZE) # fontsize of the tick labels
plt.rc('legend', fontsize=BIGGER_SIZE) # legend fontsize
plt.rc('figure', titlesize=BIGGER_SIZE) # fontsize of the figure title
#Correção de Tempo. CSFT está sempre em UTC
hours = 3
tref = '01-01-2014'
tref_file = '2014-01-01'
uba = gpd.read_file('d:\WORK\PROCAM\Dissertacao\GIS\Linha_de_Costa\LC.shp')
%matplotlib notebook
def corr_UTC(ds):
time = pd.to_datetime(ds.time.values)
time = time - timedelta(hours=hours)
ds['time'] = time
return ds
def sel_period_cal(ds):
ds = ds.where((wnd['time.year']==2014) & (ds['time.month']<6),drop=True)
return(ds)
def stickplot(T,intervalo,inicio,fim,u,v,vel,margem=5,cmap=plt.cm.YlGnBu,title='',figure_name=''):
'''
Faz a figura de stickplot
D = array/dataframe com o eixo de tempo
intervalo = intervalo de tempo (na mesma unidade do eixo de tempo) que deseja que as setas do stickplot sejam plotadas
inicio/fim = 0 e -1 se quiser plotar o período completo
= tempo em horas de início e fim que deseja cortar o período
cmap = escala de cor desejada
title = string com o título do gráfico
figure_name = string com o nome desejado para o arquivo de saída
margem = valor em horas que representa o espaço em branco no início e fim da figura para que todas as setas apareçam
'''
data_inicio = pd.Timestamp(T.year[inicio], T.month[inicio], T.day[inicio], T.hour[inicio]) - timedelta(hours=margem)
data_fim = pd.Timestamp(df.index.year[fim], df.index.month[fim], df.index.day[fim], df.index.hour[fim]) + timedelta(hours=margem)
maximo = vel[0:-1].max()
fig = plt.figure(figsize=(30,10))
ax = fig.add_subplot(1, 1, 1)
cvmap = cmap
pq = plt.quiver(T[inicio:fim:intervalo],np.zeros(len(T))[inicio:fim:intervalo],u[inicio:fim:intervalo],v[inicio:fim:intervalo],vel[inicio:fim:intervalo],
cmap = cvmap,width = 0.004,headaxislength = 3,angles = 'uv',scale_units = 'y',scale = 2.5)
plt.plot(pd.date_range(start=data_inicio, end=data_fim,freq='H'),np.zeros(len(pd.date_range(start=data_inicio, end=data_fim,freq='H'))),'--',color = 'black',linewidth = 1)
plt.grid(color='lightgray', linestyle='--', linewidth=0.8)
# Formata os eixos da figura
plt.axis([data_inicio,data_fim,-maximo-0.01,maximo+0.01])
ax.xaxis.set_major_formatter(mp_dates.DateFormatter('%d/%m %H:%M (UTC)'))
ax.yaxis.set_major_formatter(plt.NullFormatter())
fig.autofmt_xdate()
plt.title(title + '\n',fontsize = 35,weight = 'bold')
# Formata colobar
cbar = plt.colorbar(pq,cmap = cvmap)
plt.clim(0, np.ceil(maximo))
cbar.ax.set_ylabel('Intensidade (m/s)',rotation = 270,labelpad = +40)
#figname=arq.replace('.obs.p01.log','.png')
#plt.savefig(figure_name + figname, bbox_inches='tight',dpi = 200)
#plt.close('all')
def godin(wl):
a = wl.rolling(24,center=True).mean()
b = a.rolling(24,center=True).mean()
c = b.rolling(25,center=True).mean()
c = c.squeeze()
wl['Elev_Geostrof'] = c.values
wl['Elev_Astro'] = wl['water level (m)'] - c.values
return wl
def passa_baixa(sig,interval):
a = sig.rolling(interval,center=True).mean()
return a
#CFSR
wnd = xr.open_dataset('../Forcantes/CFSR/cfsr_wnd10m_2014_2020.nc')
if wnd.lon.values[0] > 180:
wnd['lon'] = wnd.lon - 360
wnd = corr_UTC(wnd)
wnd = sel_period_cal(wnd)
wnd['vel'] = ((wnd['U_GRD_L103'])**2 + (wnd['V_GRD_L103'])**2)**0.5
pres = xr.open_dataset('../Forcantes/CFSR/cfsr_prmsl_2014_2020.nc')
if pres.lon.values[0] > 180:
pres['lon'] = pres.lon - 360
pres = corr_UTC(pres)
pres = sel_period_cal(pres)
hum = xr.open_dataset('../Forcantes/CFSR/cfsr_rlhum_2014_2020.nc')
if hum.lon.values[0] > 180:
hum['lon'] = hum.lon - 360
hum = corr_UTC(hum)
hum = sel_period_cal(hum)
cloud = xr.open_dataset('../Forcantes/CFSR/cfsr_cldcovtot_2014_2020.nc')
if cloud.lon.values[0] > 180:
cloud['lon'] = cloud.lon - 360
cloud = corr_UTC(cloud)
cloud = sel_period_cal(cloud)
temp_ar = xr.open_dataset('../Forcantes/CFSR/cfsr_tmpsfc_2014_2020.nc')
if temp_ar.lon.values[0] > 180:
temp_ar['lon'] = temp_ar.lon - 360
#Convertendo de Kelvin para Celsius
temp_ar['Celsius'] = temp_ar.TMP_L1 - 273.15
temp_ar = corr_UTC(temp_ar)
temp_ar = sel_period_cal(temp_ar)
rad = xr.open_dataset('../Forcantes/CFSR/cfsr_dswsfc_2014_2020.nc')
#coordenadas
if rad.lon.values[0] > 180:
rad['lon'] = rad.lon - 360
rad = corr_UTC(rad)
rad = sel_period_cal(rad)
Nota-se que o dominio atmosférico é tem uma amplitude maior do que a área a ser considerada para a modelagem.
plt.figure(figsize = (10,10))
ax = uba.plot(color='k')
wnd.isel(time=0).vel.plot(ax=ax)
Para a validação de maré, foi utilizada apenas a estação do Porto de São Sebastião, com localização apresentada na figura a seguir:
#femar Porto de São Sebastião
fm_232 = pd.read_excel('../../dados/FEMAR/232_2014.xlsx',index_col=0,parse_dates=True)
fm_232['Amp_not_mean'] = (fm_232.Amplitude - fm_232.Amplitude.mean())/100
#femar Saco do Sombrio - Essa aqui não ficou boa
fm_231 = pd.read_excel('../../dados/FEMAR/231_2014.xlsx',index_col=0,parse_dates=True)
fm_231['Amp_not_mean'] = (fm_231.Amplitude - fm_231.Amplitude.mean())/100
#modelo
WL_SSeba = pd.read_csv('Rc_final/RC_BAROC2_FEMAR_PortoSSeba_WL.csv',index_col=0,parse_dates=True)
WL_SSeba.index = WL_SSeba.index - pd.Timedelta(hours=3)
WL_SSeba = godin(WL_SSeba)
#Validação de Nível
plt.figure(figsize=(10,10))
fm_232.Amp_not_mean.plot(color='k',linewidth=1,label='Constantes Harmônicas FEMAR Porto de São Sebastião')
WL_SSeba.Elev_Astro.plot(color='b',linewidth=0.6,linestyle=':',label='Modelo Fitrado (Godin) no FEMAR Porto de São Sebastião')
data_rec = fm_232.reindex(WL_SSeba.index)
stat = Calibration.Statistics(WL_SSeba.Elev_Astro,data_rec.Amp_not_mean)
plt.ylabel('Amplitude (m)')
plt.title('Eabs = {} Erms = {} RMAE = {}\n Skill = {} IOA = {} R2 = {}'.format(np.round(stat['Eabs'],2),np.round(stat['Erms'],2),np.round(stat['RMAE'],2),
np.round(stat['Skill'],2),np.round(stat['IOA'],2),np.round(stat['R2'],2)))
plt.xlabel('')
plt.legend()
plt.show()
#Espectro de Frequências
dt = 1 # sampling interval
Fs = 24 / dt # sampling frequency
WL_SSeba.dropna(inplace=True)
plt.figure(figsize=(10,10))
plt.magnitude_spectrum(fm_232.Amp_not_mean, Fs=Fs)
plt.title('Espectro de Frequências das Constantes Harmônicas FEMAR Porto de São Sebastião')
plt.ylabel('Amplitude (m)')
plt.xlabel('Ciclos por Dia')
plt.figure(figsize=(10,10))
plt.magnitude_spectrum(WL_SSeba['water level (m)'], Fs=Fs)
plt.title('Espectro de Frequências do Nível Total do Modelo, nas imediações da FEMAR Porto de São Sebastião')
plt.ylabel('Amplitude (m)')
plt.xlabel('Ciclos por Dia')
plt.figure(figsize=(10,10))
plt.magnitude_spectrum(WL_SSeba.Elev_Astro, Fs=Fs)
plt.title('Espectro de Frequências do Nível Astronômico do Modelo, nas imediações da FEMAR Porto de São Sebastião')
plt.ylabel('Amplitude (m)')
plt.xlabel('Ciclos por Dia')
plt.figure(figsize=(10,10))
plt.magnitude_spectrum(WL_SSeba.Elev_Geostrof, Fs=Fs)
plt.title('Espectro de Frequências do Nível Geostrófico do Modelo, nas imediações da FEMAR Porto de São Sebastião')
plt.ylabel('Amplitude (m)')
plt.xlabel('Ciclos por Dia')
Fica claro que no ponto onde foi obtido o nível total, a elevação tem forte sinal marcado em uma baixa frequência. Isto é comprovado ao passar o Filtro de Godin, pois então, o espectro de amplitude fica muito similar ao que é visto nas constantes harmônicas da FEMAR.
Conforme dito anteriormente, havia apenas uma série longa disponível para validação pontual de correntes. Segue imagem da localização do dado, coletado de 13/01/2014 - 31/05/2014
#Dados Totais do ADCP
bndo_adcp = pd.read_csv('dados/BNDO_veliger_2014.csv')
#DADOS DE CORRENTE
vel_bndo = bndo_adcp[bndo_adcp['Velocidade de Corrente [cm/s]'] != 'None']
vel_bndo.index = pd.to_datetime(vel_bndo['Data-Hora'])
vel_bndo['u'] = (vel_bndo['Velocidade de Corrente [cm/s]'].astype(float) *
np.sin(np.deg2rad(vel_bndo['Direção de Corrente [Graus]'].astype(float))))/100
vel_bndo['v'] = (vel_bndo['Velocidade de Corrente [cm/s]'].astype(float) *
np.cos(np.deg2rad(vel_bndo['Direção de Corrente [Graus]'].astype(float))))/100
#FIGURA DE COMPONENTES
ax = plt.subplots(figsize = (10,10),nrows=2,sharex=True)
vel_bndo.u.plot(ax=ax[1][0])
vel_bndo.v.plot(ax=ax[1][1])
ax[1][0].set_ylabel('Corrente (m/s)')
ax[1][1].set_ylabel('Corrente (m/s)')
ax[1][0].set_title('Componente U')
ax[1][1].set_title('Componente V')
#DADOS DE TEMPERATURA E SALINIDADE
tt = bndo_adcp[bndo_adcp['Temperatura [°C]'] != 'None']
tt.index = pd.to_datetime(tt['Data-Hora'])
ss = bndo_adcp[bndo_adcp['Salinidade [PSU]'] != 'None']
ss.index = pd.to_datetime(ss['Data-Hora'])
t0 = tt[tt['Profundidade [m]'] ==0]['Temperatura [°C]'].astype(float)
t1 = tt[tt['Profundidade [m]'] !=0]['Temperatura [°C]'].astype(float)
s0 = ss[ss['Profundidade [m]'] ==0]['Salinidade [PSU]'].astype(float)
s1 = ss[ss['Profundidade [m]'] !=0]['Salinidade [PSU]'].astype(float)
#FIGURA DE TS
ax = plt.subplots(figsize = (10,10),nrows=2,sharex=True)
t0.plot(ax=ax[1][0],label= 'Superfície')
t1.plot(ax=ax[1][0],label ='Profundidade variável')
s0.plot(ax=ax[1][1],label= 'Superfície')
s1.plot(ax=ax[1][1],label ='Profundidade variável')
ax[1][0].set_ylabel('Temperatura (°C)')
ax[1][1].set_ylabel('Salinidade (PSU)')
ax[1][0].set_title('Temperatura')
ax[1][1].set_title('Salinidade')
plt.legend()
Fica claro que no inicio da coleta há alguma inconsistência na coleta, tanto para corrente, quanto para temperatura e salinidade.
A partir do dia 11/02/2014, os dados ficam mais consistentes. Portanto para a validação feita a seguir, serão utilizados os dados apenas a partir desta data
vel_bndo = vel_bndo.loc['2014-02-11'::]
tt = tt.loc['2014-02-11'::]
ss = ss.loc['2014-02-11'::]
O modelo foi feito com 8 camadas SIGMA para tentar representar o melhor possível as estruturas termohalinas, assim sendo estratificado no verão e homogeneo no inverno.
bndo_mod_adcp = pd.read_csv('Rc_final/RC_BAROC2_BNDO_ACDP2_VEL.csv')
bndo_mod_adcp.index = pd.to_datetime(bndo_mod_adcp.iloc[:,0])
bndo_mod_adcp = bndo_mod_adcp.reindex(vel_bndo.index).interpolate()
O modelo consegue repreduzir bem o que é observado no dado, acompanhando as variações de maior período observado nos dados. Dentre as camadas, a superficial é a que melhor representa. Após esas figuras, será analisar apenas ela.
#FIGURA DE COMPONENTES
for layer in range(1,9):
ax = plt.subplots(figsize = (10,10),nrows=2,sharex=True)
print('CAMADA SIGMA {}'.format(layer))
u = bndo_mod_adcp.iloc[:,(layer*2)-1]
v = bndo_mod_adcp.iloc[:,layer*2]
vel_bndo.u.plot(ax=ax[1][0],label='Dado')
vel_bndo.v.plot(ax=ax[1][1],label='Dado')
u.plot(ax=ax[1][0],label = 'Modelo', linestyle=':')
v.plot(ax=ax[1][1],label = 'Modelo', linestyle=':')
stat_u = Calibration.Statistics(u,vel_bndo.u)
stat_v = Calibration.Statistics(v,vel_bndo.v)
ax[1][0].set_ylabel('Corrente (m/s)')
ax[1][1].set_ylabel('Corrente (m/s)')
ax[1][0].set_title('Componente U\nEabs = {} Erms = {} RMAE = {}\n Skill = {} IOA = {} R2 = {}'.format(np.round(stat_u['Eabs'],2),np.round(stat_u['Erms'],2),np.round(stat_u['RMAE'],2),
np.round(stat_u['Skill'],2),np.round(stat_u['IOA'],2),np.round(stat_u['R2'],2)))
ax[1][1].set_title('Componente V\nEabs = {} Erms = {} RMAE = {}\n Skill = {} IOA = {} R2 = {}'.format(np.round(stat_v['Eabs'],2),np.round(stat_v['Erms'],2),np.round(stat_v['RMAE'],2),
np.round(stat_v['Skill'],2),np.round(stat_v['IOA'],2),np.round(stat_v['R2'],2)))
plt.legend()
Uma leve filtragem aprimora a validação. O dado não aaparenta ter marcado a componente de maré que é visivel nos resultados do modelo.
Fica Claro nessa figura que há dois períodos distitnos registrados. Até aproximadamente dia 20/03, a maior variância no sinal da corrente é observada na Componente V. Porém, a partir desta data, este efeito passa a ser visto na Componente U.
Muito provavelmente correlacionado com o padrão de ventos.
ax = plt.subplots(figsize = (10,10),nrows=2,sharex=True)
print('CAMADA SIGMA {}'.format(layer))
u = bndo_mod_adcp.iloc[:,1]
v = bndo_mod_adcp.iloc[:,2]
u = passa_baixa(u,16)
v = passa_baixa(u,16)
vel_bndo.u.plot(ax=ax[1][0],label='Dado')
vel_bndo.v.plot(ax=ax[1][1],label='Dado')
u.plot(ax=ax[1][0],label = 'Modelo', linestyle=':')
v.plot(ax=ax[1][1],label = 'Modelo', linestyle=':')
stat_u = Calibration.Statistics(u,vel_bndo.u)
stat_v = Calibration.Statistics(v,vel_bndo.v)
ax[1][0].set_ylabel('Corrente (m/s)')
ax[1][1].set_ylabel('Corrente (m/s)')
ax[1][0].set_title('Componnente U\nEabs = {} Erms = {} RMAE = {}\n Skill = {} IOA = {} R2 = {}'.format(np.round(stat_u['Eabs'],2),np.round(stat_u['Erms'],2),np.round(stat_u['RMAE'],2),
np.round(stat_u['Skill'],2),np.round(stat_u['IOA'],2),np.round(stat_u['R2'],2)))
ax[1][1].set_title('Componnente V\Eabs = {} Erms = {} RMAE = {}\n Skill = {} IOA = {} R2 = {}'.format(np.round(stat_v['Eabs'],2),np.round(stat_v['Erms'],2),np.round(stat_v['RMAE'],2),
np.round(stat_v['Skill'],2),np.round(stat_v['IOA'],2),np.round(stat_v['R2'],2)))
plt.legend()
Na rosa de correntes, nota-se uma leve distinção da sireção e nas intensidades. Enquanto o dado aponta maior frequência de correntes para nordeste, no modelo, a direção mais frequente é a sudoeste.
No dado não é informada se foi realizada declinação dos dados, poratno é apresentada uma rosa com os dados declinados. Subtraiu-se os valores conforme calculadora da NOAA de declinação: https://www.ngdc.noaa.gov/geomag/calculators/magcalc.shtml
O dado Decliando acaba até parecendo mais do que sem declinação! Porém é um risco declinar e não saber se realmente é assim.
%matplotlib inline
HistDir(vel_bndo['Velocidade de Corrente [cm/s]'].astype(float)/100,vel_bndo['Direção de Corrente [Graus]'].astype(float),
par='corrente',Pmax=0.8,interpolado=False,logo=False)
%matplotlib inline
HistDir(vel_bndo['Velocidade de Corrente [cm/s]'].astype(float)/100,vel_bndo['Direção de Corrente [Graus]'].astype(float)-14.84,
par='corrente',Pmax=0.8,interpolado=False,logo=False)
bndo_mod_adcp['vel'] = ((bndo_mod_adcp['x component of horizontal velocity (m/s) - BNDO_ADCP_V2 - layer 1'])**2 +
(bndo_mod_adcp['y component of horizontal velocity (m/s) - BNDO_ADCP_V2 - layer 1'])**2)**0.5
bndo_mod_adcp['dir'] = np.degrees(np.arctan2(bndo_mod_adcp['x component of horizontal velocity (m/s) - BNDO_ADCP_V2 - layer 1'],
bndo_mod_adcp['y component of horizontal velocity (m/s) - BNDO_ADCP_V2 - layer 1']))
HistDir(bndo_mod_adcp.vel,bndo_mod_adcp.dir,
par='corrente',Pmax=0.8,interpolado=False,logo=False)
Nota-se que para temperatura, o modelo e o dado são mutio similares. Não tanto para o fundo, visto que é muito dificil representar a estrutura termohalina no Canal.
Para Salinidade, o modelo mostra que é mais salino que o dado. Isto provavelmente deve-se ao fato de que não coloquei inputs de água doce nem de pluviosidade.
#Lendo modelo
bndo_temp_adcp = pd.read_csv('Rc_final/RC_BAROC2_BNDO_ACDP2_TEMP.csv',encoding='latin1')
bndo_temp_adcp.index = pd.to_datetime(bndo_temp_adcp.iloc[:,0])
bndo_temp_adcp = bndo_temp_adcp.reindex(tt.index).interpolate()
bndo_sal_adcp = pd.read_csv('Rc_final/RC_BAROC2_BNDO_ACDP2_SAL.csv')
bndo_sal_adcp.index = pd.to_datetime(bndo_sal_adcp.iloc[:,0])
bndo_sal_adcp = bndo_sal_adcp.reindex(ss.index).interpolate()
#Tratando e plotando
%matplotlib notebook
t0 = tt[tt['Profundidade [m]'] ==0]['Temperatura [°C]'].astype(float)
t1 = tt[tt['Profundidade [m]'] !=0]['Temperatura [°C]'].astype(float)
s0 = ss[ss['Profundidade [m]'] ==0]['Salinidade [PSU]'].astype(float)
s1 = ss[ss['Profundidade [m]'] !=0]['Salinidade [PSU]'].astype(float)
#TEMPERATURA
ax = plt.subplots(figsize = (10,10),nrows=2,sharex=True)
t0.plot(ax=ax[1][0],label= 'Dado')
t1.plot(ax=ax[1][1],label ='Dado')
bndo_temp_adcp.iloc[:,1].plot(ax=ax[1][0],label= 'Modelo')
bndo_temp_adcp.iloc[:,8].plot(ax=ax[1][1],label= 'Modelo')
ax[1][0].set_ylabel('Temperatura (°C)')
ax[1][1].set_ylabel('Temperatura (°C)')
ax[1][0].set_title('Superfície')
ax[1][1].set_title('Fundo')
plt.legend()
#SALINIDADE
ax = plt.subplots(figsize = (10,10),nrows=2,sharex=True)
s0.plot(ax=ax[1][0],label= 'Dado')
s1.plot(ax=ax[1][1],label ='Dado')
bndo_sal_adcp.iloc[:,1].plot(ax=ax[1][0],label= 'Modelo')
bndo_sal_adcp.iloc[:,8].plot(ax=ax[1][1],label= 'Modelo')
ax[1][0].set_ylabel('Salinidade (PSU)')
ax[1][1].set_ylabel('Salinidade (PSU)')
ax[1][0].set_title('Superfície')
ax[1][1].set_title('Fundo')
plt.legend()
Aqui são algumas figuras para demonstrar a evolução temporal do modelo para outros pontos sem ser no Canal de São Sebastião.
Nas estações mais costeiras, a coluna tende a ser homogenea na maior parte do tempo. No entando, nas estações mais afastadas, observa-se a estratificação marcada quase que na maior parte do tempo. Ao final do período de modelagem é possível ver as temperaturas convergindo nestas estações, indicando o padrão esperado para a estação de inverno na plataforma.
A seguir, alguns dos pontos de monitoramento considerados. Nesta figura existem tanto pontos de coleta da CETESB, para qualidade da água costera, quanto outros dados do BNDO.
A linha em preto, indica os limites da grade.
from glob import glob
estac_temp = glob('Rc_final/*TEMP*')
for estac in estac_temp:
temp = pd.read_csv(estac,encoding='latin1')
temp.index = pd.to_datetime(temp.iloc[:,0])
plt.figure(figsize=(10,5))
temp.iloc[:,1].plot(color = 'r', label= 'Superfície - SIGMA 1')
temp.iloc[:,8].plot(color = 'cornflowerblue', label= 'Fundo - SIGMA 8 ')
plt.xlabel('')
plt.ylabel('Temperatura (°C)')
plt.title('Temperatura no Ponto {}'.format(estac.split('.')[0].split('\\')[1]))
rec = wnd.sel(lat=-23.78,lon = -45.399,method='nearest')
df = rec.to_dataframe()
No Stickplot fica evidente o mesmo padrão notado nas correntes. A partir do meio de março, ocorre uma intensificação dos ventos com mudança na direção, provavelmente indicando entradas de frente no período simulado
%matplotlib inline
cmap = plt.cm.YlGnBu
means = df.groupby(pd.Grouper(freq='3H')).mean()
stickplot(means.index,6,0,-1,means['U_GRD_L103'],means['V_GRD_L103'],df.vel.values)
means.vel.plot()
No espectro de frequência dos ventos, observa-se o padrão de baixas frequências vistos no nível do modelo, corroborando o fato de que muito provavelmente o vento é o responsável por considerável efeito na elevação e velocidade da corrente no modelo, assim como descrito na literatura
%matplotlib notebook
plt.figure(figsize=(10,5))
plt.magnitude_spectrum(rec.vel, Fs=Fs)
plt.title('Espectro de Frequências da velocidade do vento CFSR, nas imediações da FEMAR Porto de São Sebastião')
plt.ylabel('Velocidade (m/s)')
plt.xlabel('Ciclos por Dia')
Nestes gráficos fica claro como a temperatura do modelo responde a forçante atmosférica, mesmo com inputs de temperatuar e salinidade do HYCOM nas bordas. Tambéms demonstra a relaçõa com a readiação de onda curta que incide na região.
Estes efeitos mostram que o modelo realiza trocas de calor com as forçantes, mantendo um padrão esperado pelo que é forçado nele mesmo.
bndo_temp_adcp = pd.read_csv('Rc_final/RC_BAROC2_BNDO_ACDP2_TEMP.csv',encoding='latin1')
bndo_temp_adcp.index = pd.to_datetime(bndo_temp_adcp.iloc[:,0])
plt.figure(figsize=(10,5))
temp_ar.sel(lat=-23.78,lon = -45.399,method='nearest').Celsius.plot(label='Temperatura do Ar CFSR')
bndo_temp_adcp.iloc[:,1].plot(label='Modelo Superfície')
plt.title('Temperatura do Ar CFSR x Temperatura do Modelo na camada 1, nas imediações do ADCP')
plt.legend()
plt.figure(figsize=(10,5))
rad.sel(lat=-23.78,lon = -45.399,method='nearest').DSWRF_L1_Avg_1.plot()
plt.title('Radiação de Onda Curta incidente CFSR nas imediações do ADCP')